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Abstract 

This note is a physics manual for a recent numerical implementation of fct-factorization 
with running-coupling BK unintegrated gluon distributions. We also compile some results 
for Pb+Pb collisions at y/s = 2.75 TeV, such as predictions for the centrality dependence of 
the charged particle multiplicity, transverse energy, and eccentricity. The model can further 
be used to obtain initial conditions for hydrodynamic simulations of A+A collisions at the 
LHC. 

1 Introduction 

In this note we describe some of the physics underlying a recent numerical implementation of k t - 
factorization with running-coupling BK unintegrated gluon distributions for heavy-ion (and pp, 
pA) collisions at the RHIC and LHC colliders. Its main purpose is to describe the framework as 
well as some of the main assumptions and model parameters which could be tested by comparing 
to upcoming data from heavy-ion collisions at the LHC. 

Further, the model may be useful for obtaining initial conditions for hydrodynamic simulations 
of A+A collisions. Below, we also summarize a few key results. The C++ code is available for 
download at 

http : //physics .baruch. cuny . edu/node/people/adumitru/res_cgc. 

In this "manual" we do not provide an overview of the theory or applications of the Color 
Glass Condensate. 

2 The running coupling BK equation 

The Color Glass Condensate is equipped with a set of renormalization group equations, the B- 
JIMWLK equations, which describe the quantum evolution of hadron structure towards small-x 
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(see e.g. [HEiE] an d references therein). The B-JIMWLK equations are equivalent to an infinite 
set of coupled non-linear integro-differential evolution equations for the expectation values of 
the different correlators of Wilson lines averaged over the target gluon field configurations. In 
the large- A^ c limit the full B-JIMWLK hierarchy reduces to the Balitsky-Kovchegov (BK) jU [5] 
equation: a single, closed equation for the the forward scattering amplitude of a qq dipole on a 
dense target 

N{x,y,Y) = l-^(JJ{x)U\y)) Y (1) 

where the C/'s denote Wilson lines in the fundamental representation, x (y) is the transverse 
position of the quark (antiquark) and r = x — y the dipole transverse size. The average over the 
target gluon field configurations is performed at evolution rapidity Y = hx(xo/x), where x$ is the 
starting point for the evolution. The BK equation reads: 



cW(r, x) 
dlnixo/x) 



<^r_\ K(r, r 1} r 2 ) [Af(n, x) + Af(r 2 , x) - Af(r, x) - Af(n, x) Af(r 2 , x)} (2) 



where K is the evolution kernel and r 2 = r_—T\. In Eq. ([2| we have implicitly assumed translational 
invariance over scales of order the nucleon radius, i.e. mat the dipole amplitude depends only on 
the dipole transverse size r = |r| but not on the impact parameter b = (x+y)/2 of the collision^ A 
smooth variation of A/(r, x) on larger distance scales can be incorporated via the initial condition 
at the starting point x of the evolution, see below. 

The original BK equation resums small- x gluon emission to all orders at leading-logarithmic 
(LL) accuracy in a s ln(xo/a;), with a s fixed, and also contains non-linear terms that account for 
gluon-gluon self-interactions, relevant in a high density scenario. However, such limited dynamical 
input does not suffice to provide a good description of experimental data. This situation has been 
considerably improved by the recent determination of running corrections to the LL equations 
[H El IS]- Such corrections amount to just a modification of the evolution kernel in Eq. (2| with 
respect to the LL result. Though there are different possibilities of defining the running coupling 
kernel, it was shown in \X0\ that the prescription proposed by Balitsky in [8J minimizes the role of 
additional conformal corrections that arise at the same order as the running coupling, making it 
better suited for phenomenological applications. The corresponding running coupling kernel reads 
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We shall refer to Eq. (|2]) together with the evolution kernel Eq. (|3]) as the running coupling BK 
(rcBK) equation. Running coupling corrections have proven essential for promoting the BK equa- 
tion to a phenomenological tool. Indeed, the rcBK equation has been employed successfully to 
describe inclusive structure functions in e+p scattering [HI [T2], the energy and rapidity depen- 
dence of hadron multiplicities in Au+Au collisions at RHIC [13] , as well as single inclusive spectra 
in p+p and d+Au collisions at RHIC [Tlj . A detailed discussion on the numerical set up employed 
to solve the rcBK equation can be found in [10J. The present model differs from earlier predictions 
published in ref. [13] mainly due to the initial conditions described below. 

Eq. ([2]) needs to be suplemented with initial conditions. We consider two different families of 
initial conditions labeled by the parameter 7: 



Af(r,Y = 0;R) = 1 - exp 
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^^This is the main reason why we presently refrain from applying the model to compute the energy dependence 
of the multiplicity in pp collisions; accurate results require this input, see for example [B]. 
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where A = 0.241 GeV. For each of these families there are two free parameters: the value x$ 
where the evolution starts, which we fix it to Xq = 0.01 in all cases, and the initial saturation 
scale Qso(R) at the transverse coordinate R. It measures the local density of large-x sources at a 
fixed point in impact parameter space (i.e., in the transverse plane). The case 7 = 1 corresponds 
to the McLerran-Venugopalan model. For a single nucleon we take the value Q 2 = 0.2 GeV 2 
for it and we refer to it in the plots below as MV ic. However, global fits to structure functions 
in electron+proton scattering at small-x based on rcBK dynamics [TTT [T2] clearly indicate that a 
value 7 > 1 is preferred by the data. Following the results in |12^| we take the values Q 2 s0 = 0.168 
GeV 2 (for a single nucleon) and 7 = 1.119. We refer to this other set of initial conditions as 
MV 7=1 * 119 . The value of 7 determines the steepness of the unintegrated gluon distribution for 
momenta above the saturation scale k t > Q s , the MV 7=L119 initial conditions yielding a steeper 
falloff than the MV ones for all the evolution rapidities relevant at RHIC and the LHC (see fig. 1 
bottom). The main novelty in this updated version of the manual is the introduction of this new 
set of initial conditions with 7 > 1. 

As explained in more detail below, the geometry of a given A+A collision is determined by 
the fluctuations in the positions of the nucleons in the transverse plane. Each configuration 
defines a different local density in the transverse plane of each nucleus. Obviously, the smallest 
non-zero local density corresponds to the presence of a single nucleon. On the other hand, in 
A+A collisions rare fluctuations can result in collisions of a large number of nucleons at the same 
transverse position and, therefore, in a large Q s o- To account for all possible configurations we 
tabulate the solution of the rcBK equation for different values of the initial local density, i.e., for 
each value of Q s0 (R) in Eq. Q ranging from Q 2 = 0.2 (0.168) GeV 2 for MV (MV 7=1119 ) initial 
conditions (for a single nucleon) to n x Q 2 with n = 1 ... 30. The solutions are then used in the 
fct-factorization formula to calculate local gluon production at each point in the collision zone. 
Finally we perform the average over all the nucleon configurations generated by the Monte Carlo. 

To complete our discussion of the initial conditions we explain how we construct Q s q(R). 
We first generate a configuration of nucleons for each of the colliding nuclei. This consists of 
a list of random coordinates r^, i = 1 . . . A, chosen from a Woods-Saxon distribution. Multi- 
nucleon correlations are neglected except for imposing a short-distance hard core repulsion which 
enforces a minimal distance ~ 0.4 fm between any two nucleons. After this step, the longitudinal 
coordinate of any nucleon is discarded, they are projected onto the transverse plane. Factorizing 
the fluctuations of the nucleons in a nucleus from possible fluctuations of large-x "hot spots" 
within a nucleon (not accounted for at present), and finally from semi-hard gluon production 
appears to be justified by the scale hierarchy 

— < R N < Ra , (5) 

where Ra, Rn are the radii of a nucleus and of a proton, respectively. 

For a given configuration, the initial saturation momentum Q s0 (TL) at the transverse coordinate 
R is taken to be 

Q 2 (R) = iV(R)Q 2 0inucl , (6) 

where Q 2 nucl = 0.2 GeV 2 for our original uGD with MV model initial conditions while Q 2 nucl = 
0.168 GeV 2 for the new uGD with MV 7=1119 initial condition, as discussed above. AT(R) is the 
number of nucleons from the given nucleus which "overlap" the point R: 




J\T(R) = > 9 «M-|R-r,| . (7) 



2 Different parameter sets are obtained in [T2]. The values (Q 2 , 7) considered here correspond to fit h') in Table 
1, which provide a best fit to data x 2 /d.o.f = 1.104. 
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Some care must be exercised in choosing the transverse area o"o of the large-x partons of a nucleon. 
Q s0 corresponds to the density of large- x sources with x > x and should therefore be energy 
independent (recoil of the sources is neglected in the small-x approximation). We therefore take 
oo ~ 42 mb to be given by the inelastic cross-section at \/s — 200 GeV. However, <7o should not 
be confused with the energy dependent inelastic cross section a- m (s) of a nucleon which grows due 
to the emission of small-x gluons. 



3 /ct-factorization 

According to the /^-factorization formalism [15] , the number of gluons produced per unit rapidity 
at a transverse position R in A+B collisions is given by 

dN A+B^g ^ 1 da A+B->g 

dy d 2 p t d 2 R a s dy d 2 p t d 2 R ' 

where a s represents the effective interaction area and a A+B ^ 9 is the cross section for inclusive 
gluon production: 

da A+B ^ 2 1 f^tPkt f j2 , _ ,\pt + h\ n A Pt 



dy d 2 p t d 2 R 



CfpI J ~^~J 2 2 > x 2;^-&), (9) 



with x\(2) = (p t / ^/snn) exp(±y) and Cp = (N 2 — 1)/2N C . As noted before, we assume that 
the local density in each nucleus is homogenous over transverse distances of the order of the 
nucleon radius Rn- Thus, the 6-integral in Eq. ^ yields a geometric factor proportional to the 
transverse "area" of a nucleon which cancels with a similar factor implicit in a s from Eq. ([8]), 
modulo subtleties in the definition of a s . 

The unintegrated gluon distributions (ugd's) cp entering Eq. ^ are related to the dipole scat- 
tering amplitude in the adjoint representation, Mg, through a Fourier transform (for consistency 
with the notation used in Eq. (pi we make the impact parameter dependence of the ugd's explicit): 



f(k,x,b)= ° F [ d 2 re~ l]l - r V 2 r N G (r,Y = Hx /x),b). (10) 
a s (k) {2iry J 

In turn, Mg is related to the quark dipole scattering amplitude that solves the rcBK equation, A/", 
as follows: 

Af G (r,x) = 2Af(r,x) - Af 2 (r,x) . (11) 

Note that this relation entails that the saturation momentum relevant for gluon scattering is larger 
than that for quark scattering by about a factor of 2. 

Eqs. go} and Q were written originally for fixed coupling. In order to be consistent with 
our treatment of the small-x evolution, we have extended them by allowing the coupling to run 
with the momentum scale. The argument of the running coupling in Eq. d9} is chosen to be 
Q = max{|pt + kt\/2, \p t — k t \/2}, while for the definition of the ugd Eq. (\10f we take it to be 
the transverse momentum itself, k. This turns out to be important in order to reproduce the 
centrality dependence of charged particle multiplicities at RHIC, which are otherwise too flat for 
small iVp art . However, the results are not very sensitive to the particular choice of scale because 
if — > as k 2 — > due to the saturation of Af(r) at large dipole sizes r. In principle, one could 
improve on this educated ansatz by using the results of [TH] where running coupling corrections to 
inclusive gluon production have been studied. Most importantly, the x-dependence of the dipole 
scattering amplitude obtained by solving the rcBK equation encodes all the collision energy and 
rapidity dependence of the gluon production formula Eq. 
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The normalization factor K ~ 2 introduced in the /c r factorization formula (J9J) above is fixed by 
the charged particle transverse momentum distribution in p+p collisions at 7 TeV, see below. It 
lumps together higher-order corrections, sea-quark contributions, a nucleon geometry factor, and 
so on. We apply an additional "gluon multiplication factor" n g ~ 5 when computing p_i_-integrated 
yields (see below) in heavy-ion collisions but not for the transverse energy dE t /dy or for high-p 4 
hadron production in p+p from fragmenting hard gluons. 



UGD'satx=3x10" 4 . MV i.c 



3.5^ Q 2 „=0.2GeV 2 Q 2 =1.2GeV 2 Q 2 =2.4 GeV 2 




10 k t (GeV/c) 



Figure 1: Unintegrated gluon distributions for different values of the initial saturation scale evolved 
to x = 3 ■ 10~ 4 (top) for MV i.c (top). UGD's for the two different initial conditions in the single 
nucleon case: MV 7=L119 (solid) and MV (dashed) at rapidities Y = \n(xo/x) = 1.5, 3 and 6 
(bottom). 

In fig. [I] (top) we plot the ugd for three different initial MV saturation scales at x — 3 • 1CT 4 
versus transverse momentum. The ugd corresponding to a single nucleon peaks at about k t ~ 
1 GeV. The ugds for larger Q 2 s0 illustrate the shift predicted for a 6-nucleon and 12-nucleon target, 
respectively. The bottom plot shows the different k t slopes for the two initial conditions (MV 7=L119 
and MV) considered in this work in the single nucleon case. 

3.1 Observables 

Eq. ([9]) is the starting point for all observables shown below. In particular, the charged particle 
multiplicity and the transverse energy can be obtained by integrating over the transverse plane 
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and p t , 



dN ch 2 f 2 f 2 dN A+B ^ 

K g d R d p t „ 2 (12) 



cfy 3 y 7 J rL dyd 2 p t d 2 R 

d R — D h (z = -) -j—j^—j^- (13) 



d?p t dy J J z 2 k t dy d 2 k t d 2 R 

1 d 2 R d 2 p t p t -———. (14) 



dy J J dy d 2 p t d 2 R 

Note that a low-p^ cutoff is not required since the integration over fc t in d9| extends only up to 
Pt- The saturation of the gluon distribution functions guarantees that the dominant scale in the 
transverse momentum integrations is the saturation momentum. For the single-inclusive hadron 



Pt distributions (13) we use the KKP gluon — > charged hadron LO fragmentation function [T7] 
with the scale Q^= k 2 ; the integral over the hadron momentum fraction is restricted to z > 0.05 
to avoid a violation of the momentum sum rule. The "gluon multiplication factor" is fixed to 
Kg = 5 (energy and centrality independent) in order to reproduce the measured charged hadron 
multiplicity in heavy-ion collisions at RHIC and LHC energies, see below. The upper limit in 
the integrals over the gluon transverse momentum in dN^/dy and dE t /dy has been taken as 
£>max = 12 GeV; if the integrals are extended further then a slight adjustment of the normalization 
factors K g and K may be required. 
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Figure 2: Centrality dependence of the charged particle multiplicity at midrapidity for Au+Au 
collisions at = 200 GeV. PHOBOS data: ref. [IE]. 

In order to compare our results for initial gluon production to the final state distributions of 
detected particles one has to translate the rapidity distributions into pseudo-rapidity distributions 
through the y — >■ rj Jacobian, 



dN, 



ch 



cosily 



dN, 



ch 



cLE^ 
drj 



dy 



cosh 2 f] + m 2 / P 2 
cosh 7] dE t 

cosh 2 rj + m 2 / P 2 y 



(15) 
(16) 
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Figure 3: Centrality dependence of the charged particle multiplicity at midrapidity for Pb+Pb 
collisions at y/s = 2.76 TeV. Alice data from ref. 



with y = | In (y cosh 2 r\ + m 2 / P 2 + sinh?7)/(y cosh 2 r\ + m 2 /P 2 — sinhr/). We assume that in 

this Jacobian m = 350 MeV and P = 0.13 GeV + 0.32 GeV^/s/1 TeV°' U5 . Note that such 
transformation is not needed (it is trivial) if one is interested in initial (massless) gluon production 
to initialize a hydrodynamic simulation. 



4 Results 

In this section we present some results obtained with the new MV 7=1119 uGD and compare them 
to the original uGD obtained by rcBK evolution with MV model initial conditions. 

In Fig. [2] we show the centrality dependence of dN/dr] at full RHIC energy. This energy probes 
mainly the initial condition for small-x evolution. Both sets perform reasonably well and we do 
not find it possible to distinguish. 

Fig. [3 shows the centrality dependence of dN/dr] at 77 = for Pb+Pb collisions at y/s = 
2.76 TeVfl Once again both uGD sets agree with the data to within a few percent. It is worth 
noting that our present results for the multiplicity in central collisions are higher than the one 
predicted in ref. [IB], which also relied on the rcBK equation plus fct-factorization. The main 
difference between these two works relates to the treatment of the geometry of the collision: while 
in [13] the full nucleus was characterized by a single average saturation scale and then evolved 
to higher energies, the Monte Carlo approach presented here allows for different local densities at 
every point in the transverse plane, each of them evolved locally to higher energies. The average 
over different configurations is performed after the evolution, and not before, as implicitly done in 
[T3] . Thus we interpret these two different results as an indication that the average over nuclear 
geometry does not commute with the evolution. 

In Fig. [4] we show the transverse momentum distribution of charged particles for p+p colli- 
sions at ^/s — 7 TeV. We recall that presently we assume uniform and homogeneous, "disc-like" 
nucleons. However, we do not attempt to compute the inelastic proton-proton cross section and 
high-p 4 hadron spectra may not be very sensitive to the proton impact parameter profile. 



3 The Alice papers summarize other recent predictions for dN c h/drj for Pb+Pb collisions. 
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Figure 4: Transverse momentum distribution of charged particles at r] = for p+p collisions at 
y/s = 7 TeV. CMS data from ref. [20]. 



For the range of p± shown in the figure, particle production probes LC momentum fractions 
well below our assumed starting point of Xq = 0.01. The uGD derived from MV model initial 
conditions is clearly too "hard" and predicts an incorrect slope. The new uGD obtained from the 
MV 7=1119 initial condition corrects this deficiency and provides a good description of the CMS 
data in the small-x, semi-hard regime (see, also, ref. j6j). This illustrates the power of LHC to 
constrain small- x physics. Also, we have used this observable to fix the genuine "K-factor" to 
K = 2 (MV 7=L119 i.e.) or K = 1.5 (MV i.e.), respectively. 
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Figure 5: Transverse energy at rj = versus centrality for Pb+Pb collisions at -^/s 
(left) and for Au+Au collisions at y/s = 200 GeV (right). 



2.76 TeV 



The transverse energy in heavy-ion collisions versus centrality is shown in Fig. [5j The new 
MV 7=1119 uGD predicts a lower transverse energy than the one derived from MV model initial 
conditions, by about 1 GeV per participant at LHC energies^} It should be kept in mind that the 
measured transverse energy in the final state is only a lower bound on the initial dE±/drj since soft 
final-state interactions such as collective longitudinal flow do not conserve the transverse energy 

4 Note that we deviate from vl of this manuscript in that the "gluon multiplication factor" k 9 has been dropped 
from the normalization for Ej_. 
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Figure 6: Charged particle pseudo-rapidity distributions for UGD Set 1 (MV model initial con- 
ditions) in pp collisions at 7 TeV (normalization adjustement -10%), 2.36 TeV (normalization 
adjustement -4%), and p+Pb collisions at 4.4 TeV (standard normalization as for AA collisions 
from Figs. |2l II. 



In Figj6]we show the rapidity distributions in pp collisions at two LHC energies and a prediction 
for minimum bias p+Pb collisions at 4.4 TeV. We restrict these distributions to \r]\ < 2.5 since we 
have not accounted for the contribution of scattered valence partons. While there is perhaps an 
overall normalization uncertainty of up to ~ 10%, we note the large forward-backward asymmetry 
for p+Pb collisiong^it rapidities far from the valence partons predicted by the rcBK-UGD. 

We leave a more detailed study of the sensitivity of our results to the parameters of the model 
and a more detailed comparison with RHIC and LHC data for a future publication. 
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